home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / p_correlate.pro < prev    next >
Text File  |  1997-07-08  |  4KB  |  112 lines
  1. ;$Id: p_correlate.pro,v 1.6 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       P_CORRELATE
  8. ;
  9. ; PURPOSE:
  10. ;       This function computes the partial correlation coefficient of a
  11. ;       dependent variable and one particular independent variable when
  12. ;       the effects of all other variables involved are removed.
  13. ;
  14. ; CATEGORY:
  15. ;       Statistics.
  16. ;
  17. ; CALLING SEQUENCE: 
  18. ;       Result = P_correlate(X, Y, C)
  19. ;
  20. ; INPUTS:
  21. ;       X:    An n-element vector of type integer, float or double that
  22. ;             specifies the independent variable data.
  23. ;       
  24. ;       Y:    An n-element vector of type integer, float or double that
  25. ;             specifies the dependent variable data.
  26. ;
  27. ;       C:    An array of type integer, float or double that specifies the 
  28. ;             independent variable data whose effects are to be removed. 
  29. ;             The columns of this two dimensional array correspond to the 
  30. ;             n-element vectors of independent variable data.
  31. ;
  32. ; KEYWORD PARAMETERS:
  33. ;    DOUBLE:  If set to a non-zero value, computations are done in
  34. ;             double precision arithmetic.
  35. ;
  36. ; EXAMPLES:
  37. ;       Define the data vectors.
  38. ;         x0 = [64, 71, 53, 67, 55, 58, 77, 57, 56, 51, 76, 68]
  39. ;         x1 = [57, 59, 49, 62, 51, 50, 55, 48, 52, 42, 61, 57]
  40. ;         x2 = [ 8, 10,  6, 11,  8,  7, 10,  9, 10,  6, 12,  9]
  41. ;
  42. ;       Compute the partial correlation of x0 and x1 with the effects of
  43. ;       x2 removed. The result should be 0.533469
  44. ;         result = p_correlate(x0, x1, reform(x2, 1, n_elements(x2)))
  45. ;
  46. ;       Compute the partial correlation of x0 and x2 with the effects of 
  47. ;       x1 removed. The result should be 0.334572
  48. ;         result = p_correlate(x0, x2, reform(x1, 1, n_elements(x1)))
  49. ;
  50. ;       Compute the partial correlation of x1 and x2 with the effects of 
  51. ;       x0 removed. The result should be 0.457907
  52. ;         result = p_correlate(x1, x2, reform(x0, 1, n_elements(x0)))
  53. ;
  54. ; REFERENCE:
  55. ;       APPLIED STATISTICS (third edition)
  56. ;       J. Neter, W. Wasserman, G.A. Whitmore
  57. ;       ISBN 0-205-10328-6
  58. ;
  59. ; MODIFICATION HISTORY:
  60. ;       Modified by:  GGS, RSI, July 1994
  61. ;                     Minor changes to code. New documentation header.
  62. ;       Modified by:  GGS, RSI, August 1996
  63. ;                     Added DOUBLE keyword. 
  64. ;                     Modified keyword checking and use of double precision.
  65. ;-
  66.  
  67. FUNCTION P_Correlate, X, Y, C, Double = Double
  68.  
  69.   ON_ERROR, 2  ;Return to caller if an error occurs.
  70.  
  71.   Sx = SIZE(x)  &  Sy = SIZE(y)  &  Sc = SIZE(c)
  72.  
  73.   if Sx[Sx[0]+2] ne Sy[Sy[0]+2] then MESSAGE, $
  74.     "X and Y must have the same number of elements."
  75.  
  76.   if Sc[0] ne 2 then MESSAGE, $
  77.     "C parameter must be a two-dimensional array."
  78.  
  79.   ;Check row dimension of C.
  80.   if Sx[Sx[0]+2] ne Sc[Sc[0]] then MESSAGE, $
  81.     "Incompatible arrays."
  82.  
  83.   if N_ELEMENTS(Double) eq 0 then $
  84.     Double = (Sx[Sx[0]+1] eq 5) or (Sy[Sy[0]+1] eq 5) or (Sc[Sc[0]+1] eq 5)
  85.  
  86.   if Sc[1] eq 1 then begin
  87.     p = [CORRELATE(X, Y, Double = Double), $
  88.          CORRELATE(X, C, Double = Double), $
  89.          CORRELATE(Y, C, Double = Double)]
  90.     if (p[1] ne 1 and p[2] ne 1) then $
  91.       RETURN, (p[0] - p[1] * p[2])/SQRT((1 - p[1]^2) * (1 - p[2]^2)) $
  92.       else RETURN, 0 * p
  93.   endif else begin
  94.     ;Vector of weights.
  95.     if Double eq 0 then Wts = REPLICATE(1.0, Sc[2]) $
  96.     else Wts = REPLICATE(1.0d, Sc[2])
  97.     cf = REGRESS(C, Y, Wts, yFit, a0, s, f, r, p0, /RELATIVE)
  98.     cf = REGRESS([C, TRANSPOSE(X)], Y, Wts, yFit, a0, s, f, r, p1, /RELATIVE)
  99.     if Double eq 0 then begin
  100.       p0 = FLOAT(p0)
  101.       p1 = FLOAT(p1)
  102.     endif
  103.     p0 = 1 - p0^2
  104.     p1 = 1 - p1^2 
  105.     if p0 eq 0 then $
  106.       RETURN, 0 * p $
  107.       else RETURN, SQRT((p0 - p1)/p0) 
  108.   endelse
  109.  
  110. END
  111.  
  112.